This document serves as a supplementary material for my thesis demonstrating some of the scripting done to achieve the presented results. It should allow anyone to reproduce the data and evaluate the results for themselves.
## Warning: package 'flowCore' was built under R version 4.1.1
## Warning: package 'rjson' was built under R version 4.1.2
## Warning: package 'flowWorkspace' was built under R version 4.1.1
## As part of improvements to flowWorkspace, some behavior of
## GatingSet objects has changed. For details, please read the section
## titled "The cytoframe and cytoset classes" in the package vignette:
##
## vignette("flowWorkspace-Introduction", "flowWorkspace")
## Warning: package 'FNN' was built under R version 4.1.3
## Warning: package 'ggplot2' was built under R version 4.1.3
## Warning: package 'EmbedSOM' was built under R version 4.1.3
## Warning: package 'gridExtra' was built under R version 4.1.2
## Warning: package 'cowplot' was built under R version 4.1.2
## Warning: package 'ggpubr' was built under R version 4.1.3
##
## Attaching package: 'ggpubr'
## The following object is masked from 'package:cowplot':
##
## get_legend
## Warning: package 'scales' was built under R version 4.1.3
Set working directory to the location of this notebook.
try({
dr = getSrcDirectory()[1]
})
## Error in getSrcFilename(x, full.names = TRUE, unique = unique) :
## argument "x" is missing, with no default
#when sourcing the file
try({
dr = dirname(rstudioapi::getActiveDocumentContext()$path)
})
#for running code directly in rstudio
setwd(dr) #set for usage as a script/running lines
knitr::opts_knit$set(root.dir = dr) #set for R markdown chunk execution
spctr = fromJSON(file = "panel1_lymph_subtracted_fixed.json")
spectra = matrix(ncol = 64, nrow = length(spctr))
spec_names = c()
for (i in 1:length(spctr)) {
spectra[i, ] = spctr[[i]]$spectrum$mS
spec_names = c(spec_names, spctr[[i]]$antigen)
}
colnames(spectra) = spctr[[1]]$spectrum$channels
rm(spctr)
rownames(spectra) = spec_names
nspectra = nrow(spectra)
ndetects = ncol(spectra)
Unpack the tree-like structure of the phenotype table using previously defined recursive function.
pheno = read.csv('phenotypes_noAFonly.csv',
sep = ',',
fileEncoding = "UTF-8-BOM")
todel = c()
pheno = recur_unpack(c('base'), pheno, todel)
#write.csv(pheno,'pheno_unpacked.csv') #if you are interested in unpacked phenotypes uncomment this
Change nthreads parameter appropriately up to the number of threads available to your machine.
stronk = c(2, 3, 9, 14, 17, 21)
res = generate_artificial_cytometry_data(spectra, pheno, 5000, stronk)
received = res[[1]]
gt = res[[2]]
nr = res[[3]]
len_umx = length(gt)
len_rec = length(received)
start.time = Sys.time()
ngd = nougadmt(
received,
spectra = spectra,
snw = 5,
spw = 1,
nw = 200,
start = 2,
iters = 500,
alpha = 0.01,
nthreads = detectCores()
)$unmixed
stop.time = Sys.time()
cat('Multithreaded nougad: ')
## Multithreaded nougad:
stop.time - start.time
## Time difference of 0.498981 secs
If you want to compare performance with the singlethreaded implementation.
start.time = Sys.time()
ngd = nougad(
received,
spectra = spectra,
snw = 5,
spw = 1,
nw = 200,
start = 2,
iters = 500,
alpha = 0.01
)$unmixed
stop.time = Sys.time()
cat('Singlehreaded nougad: ')
## Singlehreaded nougad:
stop.time - start.time
## Time difference of 4.479315 secs
(Speed up from going multithreaded is almost 20x on my Ryzen 9 5900X CPU with 12cores/24threads and 64GB of RAM. For 100000 cells this means going from ~112s to ~5.7s which is significant.)
ols = t(lm(t(received) ~ t(spectra) + 0)$coefficients)
wols = ols
for (i in 1:nr) {
ws = norm(pmax(0, received[i, ] / sum(received[i, ])))
wols[i, ] = t(lm(t(received[i, , drop = F]) ~ t(spectra) + 0, weights =
ws)$coefficients)
}
ngd_t = trans(ngd)
ols_t = trans(ols)
wols_t = trans(wols)
gt_t = trans(gt)
ngd_diff = mse(gt_t, ngd_t)
wols_diff = mse(gt_t, wols_t)
ols_diff = mse(gt_t, ols_t)
cat(paste0('MSE for NGD: ', ngd_diff , '\n'))
## MSE for NGD: 0.762993576176325
cat(paste0('MSE for WOLS: ', wols_diff , '\n'))
## MSE for WOLS: 1.88047084606996
cat(paste0('MSE for OLS: ', ols_diff , '\n'))
## MSE for OLS: 1.24559027180041
ngd_err = (ngd_t - gt_t) ^ 2
ols_err = (ols_t - gt_t) ^ 2
wols_err = (wols_t - gt_t) ^ 2
x = 1:len_umx
par(
mfrow = c(3, 1),
tcl = 0.5,
mgp = c(0, -1.4, 0),
mar = c(0, 0, 1.5, 0)
)
plot(
x,
ngd_err,
type = 'l',
col = 'red',
pch = 16,
cex = .1,
ylim = c(0, 50),
main = paste0('NOUGAD MSE=', round(ngd_diff, 4), ' SD_mse=', round(sd(ngd_err), 4)),
xaxt = 'n'
)
plot(
x,
wols_err,
type = 'l',
col = 'red',
pch = 16,
cex = .1,
ylim = c(0, 50),
main = paste0('WOLS MSE=', round(wols_diff, 4), ' SD_mse=', round(sd(wols_err), 4)),
xaxt = 'n'
)
plot(
x,
ols_err,
type = 'l',
col = 'red',
pch = 16,
cex = .1,
ylim = c(0, 50),
main = paste0('OLS MSE=', round(ols_diff, 4), ' SD_mse=', round(sd(ols_err), 4)),
xaxt = 'n'
)
res_1k = generate_artificial_cytometry_data(spectra, pheno, 1000, stronk)
received_1k = res_1k[[1]]
gt_1k = res_1k[[2]]
nr_1k = res_1k[[3]]
ngd_1k = nougadmt(
received_1k,
spectra = spectra,
snw = 5,
spw = 1,
nw = 200,
start = 2,
iters = 500,
alpha = 0.01,
nthreads=detectCores()
)$unmixed
ols_1k = t(lm(t(received_1k) ~ t(spectra) + 0)$coefficients)
wols_1k = ols_1k
for (i in 1:nr_1k) {
ws = norm(pmax(0, received_1k[i,] / sum(received_1k[i,])))
wols_1k[i,] = t(lm(t(received_1k[i, , drop = F]) ~ t(spectra) + 0, weights =
ws)$coefficients)
}
ngd_t_1k = trans(ngd_1k)
ols_t_1k = trans(ols_1k)
wols_t_1k = trans(wols_1k)
gt_t_1k = trans(gt_1k)
ngd_err_1k = (ngd_t_1k - gt_t_1k) ^ 2
ols_err_1k = (ols_t_1k - gt_t_1k) ^ 2
wols_err_1k = (wols_t_1k - gt_t_1k) ^ 2
total_brightness_1k = rowSums(received_1k)
ngd_mcr_1k = rowSums(ngd_err_1k) / nspectra
ols_mcr_1k = rowSums(ols_err_1k) / nspectra
wols_mcr_1k = rowSums(wols_err_1k) / nspectra
rnms_1k = rownames(received_1k)
rnms_1k = sapply(strsplit(rnms_1k, split = '.', fixed = TRUE), function(x)
(x[1]))
uq = unique(rnms_1k)
occurs = Vectorize(grep, 'pattern')(paste0('^', uq, '$'), rnms_1k)
cnts = as.vector(sapply(occurs, length))
coords = c(0)
for (i in 1:length(cnts))
{
coords[i + 1] = sum(cnts[1:i])
}
coords = coords[-length(coords)]
end_point = 0.5 + length(uq) + length(uq) - 1 #this is the line which does the trick (together with barplot "space = 1" parameter)
x = 1:nr_1k
par(
mfrow = c(3, 1),
tcl = 0.5,
mgp = c(0, -1.4, 0),
mar = c(0, 0, 1.5, 0)
)
plot(
x,
ngd_mcr_1k,
pch = 16,
cex = .1,
ylim = c(0, 15),
main = 'NOUGAD',
type = 'l',
col = 'red',
xaxt = 'n',
yaxt = 'n'
)
lines(x, norm(total_brightness_1k) * 15 , col = alpha('green', 0.5))
legend(
round(nr_1k - 1 / 7 * nr_1k),
15,
legend = c("Squared error", "Cell luminance"),
col = c("red", "green"),
lty = 1:2
)
axis(2,
at = seq(0, 15, 2),
labels = T,
gap.axis = 0)
plot(
x,
wols_mcr_1k,
pch = 16,
cex = .1,
ylim = c(0, 15),
main = 'WOLS',
type = 'l',
col = 'red',
xaxt = 'n',
xlab = '',
yaxt = 'n'
)
lines(x, norm(total_brightness_1k) * 15 , col = alpha('green', 0.5))
text(
coords,
par("usr")[3] + 15,
srt = 90,
adj = 1,
xpd = TRUE,
labels = uq,
cex = 0.9
)
axis(1,
gap.axis = 1,
padj = 0.5,
xaxt = 'n')
axis(2,
at = seq(0, 15, 2),
labels = T,
gap.axis = 0)
plot(
x,
ols_mcr_1k,
pch = 16,
cex = .1,
ylim = c(0, 15),
main = 'OLS',
type = 'l',
col = 'red',
xaxt = 'n',
yaxt = 'n'
)
lines(x, norm(total_brightness_1k) * 15 , col = alpha('green', 0.5))
axis(2,
at = seq(0, 15, 2),
labels = T,
gap.axis = 0)
Correlation of mean per cell error with total cell brightness:
ngd_err = (ngd_t - gt_t) ^ 2
ols_err = (ols_t - gt_t) ^ 2
wols_err = (wols_t - gt_t) ^ 2
total_brightness = rowSums(received)
ngd_mcr = rowSums(ngd_err) / nspectra
ols_mcr = rowSums(ols_err) / nspectra
wols_mcr = rowSums(wols_err) / nspectra
df_ngd = data.frame(
cell_brightness = norm(total_brightness),
mean_cell_error = norm(ngd_mcr)
)
df_ngd = df_ngd[order(df_ngd$mean_cell_error),]
df_ols = data.frame(
cell_brightness = norm(total_brightness),
mean_cell_error = norm(ols_mcr)
)
df_ols = df_ols[order(df_ols$mean_cell_error),]
df_wols = data.frame(
cell_brightness = norm(total_brightness),
mean_cell_error = norm(wols_mcr)
)
df_wols = df_wols[order(df_wols$mean_cell_error),]
cn = cor.test(df_ngd$cell_brightness, df_ngd$mean_cell_error)
co = cor.test(df_ols$cell_brightness, df_ols$mean_cell_error)
cw = cor.test(df_wols$cell_brightness, df_wols$mean_cell_error)
par(mar = c(0, 0, 3, 0), mfrow = c(1, 3))
plot(
df_ngd$cell_brightness,
df_ngd$mean_cell_error,
pch = 16,
cex = .1,
main = paste0(
'NOUGAD r=',
round(cn$estimate, 3),
' 95CI=(',
round(cn$conf.int[1], 2),
';',
round(cn$conf.int[2], 2),
')\n p',
format.pval(cn$p.value)
),
xaxt = 'n',
yaxt = 'n'
)
abline(0, 1, col = 'red')
plot(
df_wols$cell_brightness,
df_wols$mean_cell_error,
pch = 16,
cex = .1,
main = paste0(
'WOLS r=',
round(cw$estimate, 3),
' 95CI=(',
round(cw$conf.int[1], 2),
';',
round(cw$conf.int[2], 2),
')\n p',
format.pval(cw$p.value)
),
xaxt = 'n',
yaxt = 'n'
)
abline(0, 1, col = 'red')
plot(
df_ols$cell_brightness,
df_ols$mean_cell_error,
pch = 16,
cex = .1,
main = paste0(
'OLS r=',
round(co$estimate, 3),
' 95CI=(',
round(co$conf.int[1], 2),
';',
round(co$conf.int[2], 2),
')\n p',
format.pval(co$p.value)
),
xaxt = 'n',
yaxt = 'n'
)
abline(0, 1, col = 'red')
spec_int = norm(rowSums(spectra) ^ 2) #Total energy of the given spectrum
mean_ag_expr = norm(colSums(gt) / nspectra)#Mean expression of that spectrum in the data (ground truth)
colvec <- colorRampPalette(c("green", "yellow", "red"))
cvc = c(colvec(10), rep('red', 10))
for (i in 1:nspectra) {
par(mar = c(0, 2, 2, 0), mfrow = c(1, 3))
ols_vec = ols_err[, i]
wols_vec = wols_err[, i]
ngd_vec = ngd_err[, i]
xmx = max(c(max(ols_vec), max(wols_vec), max(ngd_vec)))
h_ols =
hist(ols_vec,
plot = FALSE,
breaks = seq(0, xmx, length.out = 20))
h_ols$counts = log(h_ols$counts)
h_ols$counts[h_ols$counts == -Inf] = 0
h_wols =
hist(wols_vec,
plot = FALSE,
breaks = seq(0, xmx, length.out = 20))
h_wols$counts = log(h_wols$counts)
h_wols$counts[h_wols$counts == -Inf] = 0
h_ngd =
hist(ngd_vec,
plot = FALSE,
breaks = seq(0, xmx, length.out = 20))
h_ngd$counts = log(h_ngd$counts)
h_ngd$counts[h_ngd$counts == -Inf] = 0
ymx = max(c(max(h_ols$counts), max(h_wols$counts), max(h_ngd$counts)))
plot(
h_ngd,
col = cvc,
main = paste0(
spec_names[i],
' NGD\n E=',
round(spec_int[i], 3),
' M(Ex)=',
round(mean_ag_expr[i], 3)
),
xaxt = 'n',
xlim = c(0, xmx),
ylim = c(0, ymx)
)
plot(
h_wols,
col = cvc,
main = paste0(
spec_names[i],
' WOLS\n E=',
round(spec_int[i], 3),
' M(Ex)=',
round(mean_ag_expr[i], 3)
),
xaxt = 'n',
xlim = c(0, xmx),
ylim = c(0, ymx),
yaxt = 'n'
)
plot(
h_ols,
col = cvc,
main = paste0(
spec_names[i],
' OLS\n E=',
round(spec_int[i], 3),
' M(Ex)r=',
round(mean_ag_expr[i], 3)
),
xaxt = 'n',
xlim = c(0, xmx),
ylim = c(0, ymx),
yaxt = 'n'
)
}
Correlation of mean spectrum error with the spectrum energy.
ngd_msr = colSums(ngd_err) / nr
ols_msr = colSums(ols_err) / nr
wols_msr = colSums(wols_err) / nr
df_ngd = data.frame(spectrum_energy = norm(spec_int),
mean_spectrum_error = norm(ngd_msr))
df_ngd = df_ngd[order(df_ngd$mean_spectrum_error),]
df_ols = data.frame(spectrum_energy = norm(spec_int),
mean_spectrum_error = norm(ols_msr))
df_ols = df_ols[order(df_ols$mean_spectrum_error),]
df_wols = data.frame(spectrum_energy = norm(spec_int),
mean_spectrum_error = norm(wols_msr))
df_wols = df_wols[order(df_wols$mean_spectrum_error),]
cn = cor.test(ngd_msr, spec_int)
co = cor.test(ols_msr, spec_int)
cw = cor.test(wols_msr, spec_int)
par(mar = c(0, 0, 3, 0), mfrow = c(1, 3))
plot(
df_ngd$spectrum_energy,
df_ngd$mean_spectrum_error,
pch = 16,
cex = .1,
main = paste0(
'NOUGAD r=',
round(cn$estimate, 3),
' 95CI=(',
round(cn$conf.int[1], 2),
';',
round(cn$conf.int[2], 2),
')\n p=',
format.pval(cn$p.value)
),
xaxt = 'n',
yaxt = 'n'
)
abline(0, 1, col = 'red')
plot(
df_wols$spectrum_energy,
df_wols$mean_spectrum_error,
pch = 16,
cex = .1,
main = paste0(
'WOLS r=',
round(cw$estimate, 3),
' 95CI=(',
round(cw$conf.int[1], 2),
';',
round(cw$conf.int[2], 2),
')\n p=',
format.pval(cw$p.value)
),
xaxt = 'n',
yaxt = 'n'
)
abline(0, 1, col = 'red')
plot(
df_ols$spectrum_energy,
df_ols$mean_spectrum_error,
pch = 16,
cex = .1,
main = paste0(
'OLS r=',
round(co$estimate, 3),
' 95CI=(',
round(co$conf.int[1], 2),
';',
round(co$conf.int[2], 2),
')\n p=',
format.pval(co$p.value)
),
xaxt = 'n',
yaxt = 'n'
)
abline(0, 1, col = 'red')
Using lines.
colnames(gt_t_1k) = spec_names
colnames(ngd_t_1k) = spec_names
colnames(ols_t_1k) = spec_names
colnames(wols_t_1k) = spec_names
ngd_t_1k = as.data.frame(ngd_t_1k)
gt_t_1k = as.data.frame(gt_t_1k)
ols_t_1k = as.data.frame(ols_t_1k)
wols_t_1k = as.data.frame(wols_t_1k)
marker_cords = which(mean_ag_expr > 0.3)
markers = spec_names[marker_cords]
markers = markers[markers != 'Autofluorescence']
combos = combn(markers, 2)
for (i in 1:ncol(combos)) {
df_ngd = data.frame(
X1 = gt_t_1k[combos[1, i]][, 1],
X2 = gt_t_1k[combos[2, i]][, 1],
X1.1 = ngd_t_1k[combos[1, i]][, 1],
X2.1 = ngd_t_1k[combos[2, i]][, 1]
)
df_ols = data.frame(
X1 = gt_t_1k[combos[1, i]][, 1],
X2 = gt_t_1k[combos[2, i]][, 1],
X1.1 = ols_t_1k[combos[1, i]][, 1],
X2.1 = ols_t_1k[combos[2, i]][, 1]
)
df_wols = data.frame(
X1 = gt_t_1k[combos[1, i]][, 1],
X2 = gt_t_1k[combos[2, i]][, 1],
X1.1 = wols_t_1k[combos[1, i]][, 1],
X2.1 = wols_t_1k[combos[2, i]][, 1]
)
p1 = ggplot() + geom_segment(
data = df_ngd,
aes(
x = X1,
y = X2,
xend = X1.1,
yend = X2.1
),
colour = "red",
alpha = 0.5
) + geom_point(aes(x = df_ngd$X1, y = df_ngd$X2),
size = 0.5,
alpha = 0.3) + theme_cowplot() + ylab(combos[2, i]) + xlab('') + ggtitle('NGD')
p2 = ggplot() + geom_segment(
data = df_wols,
aes(
x = X1,
y = X2,
xend = X1.1,
yend = X2.1
),
colour = "red",
alpha = 0.5
) + geom_point(aes(x = df_ngd$X1, y = df_ngd$X2),
size = 0.5,
alpha = 0.3) + theme_cowplot() + ylab('') + xlab(combos[1, i]) + ggtitle('WOLS')
p3 = ggplot() + geom_segment(
data = df_ols,
aes(
x = X1,
y = X2,
xend = X1.1,
yend = X2.1
),
colour = "red",
alpha = 0.5
) + geom_point(aes(x = df_ngd$X1, y = df_ngd$X2),
size = 0.5,
alpha = 0.3) + theme_cowplot() + ylab('') + xlab('') + ggtitle('OLS')
print(ggarrange(p1, p2, p3, ncol = 3))
}
Color cells by squared error magnitude.
colnames(gt_t) = spec_names
colnames(ngd_t) = spec_names
colnames(ols_t) = spec_names
colnames(wols_t) = spec_names
ngd_t = as.data.frame(ngd_t)
gt_t = as.data.frame(gt_t)
ols_t = as.data.frame(ols_t)
wols_t = as.data.frame(wols_t)
ngd_err = as.data.frame(ngd_err)
colnames(ngd_err) = spec_names
ols_err = as.data.frame(ols_err)
colnames(ols_err) = spec_names
wols_err = as.data.frame(wols_err)
colnames(wols_err) = spec_names
for (i in 1:ncol(combos)) {
aes_engd = rowSums(ngd_err[, c(combos[1, i], combos[2, i])]) / 2
aes_eols = rowSums(ols_err[, c(combos[1, i], combos[2, i])]) / 2
aes_ewols = rowSums(wols_err[, c(combos[1, i], combos[2, i])]) / 2
mx = max(c(max(aes_engd), max(aes_eols), max(aes_ewols)))
x_gt = gt_t[combos[1, i]][, 1]
y_gt = gt_t[combos[2, i]][, 1]
x_ngd = ngd_t[combos[1, i]][, 1]
y_ngd = ngd_t[combos[2, i]][, 1]
x_ols = ols_t[combos[1, i]][, 1]
y_ols = ols_t[combos[2, i]][, 1]
x_wols = wols_t[combos[1, i]][, 1]
y_wols = wols_t[combos[2, i]][, 1]
mx_x = max(c(max(x_ols), max(x_wols), max(x_ngd), max(x_gt)))
mx_y = max(c(max(y_ols), max(y_wols), max(y_ngd), max(y_gt)))
mn_x = min(c(min(x_ols), min(x_wols), min(x_ngd), min(x_gt)))
mn_y = min(c(min(y_ols), min(y_wols), min(y_ngd), min(y_gt)))
Error = rep(0, nr)
p1 = qplot() + geom_point(aes(x_gt, y_gt, colour = Error),
size = 1,
alpha = 0.5) + scale_colour_gradientn(
limits = c(0, mx),
trans = 'pseudo_log',
colours = c("blue", "yellow", "red"),
labels = c(' 0', ' 2', ' 5', '10 ', '20 ', '40 ', '60 '),
breaks = c(0, 2, 5, 10, 20, 40, 60)
) + ggtitle('Ground Truth') + theme_cowplot() +
ylab('') + xlab('') + scale_x_continuous(limits = c(mn_x, mx_x)) + scale_y_continuous(limits = c(mn_y, mx_y)) +
theme_void()
p2 = qplot() + geom_point(aes(x_ngd, y_ngd, colour = aes_engd),
size = 1,
alpha = 0.5) + scale_colour_gradientn(
limits = c(0, mx),
trans = 'pseudo_log',
colours = c("blue", "yellow", "red"),
labels = c(' 0', ' 2', ' 5', '10 ', '20 ', '40 ', '60 '),
breaks = c(0, 2, 5, 10, 20, 40, 60)
) + ggtitle('NGD') + theme_cowplot() +
ylab('') + xlab('') + scale_x_continuous(limits = c(mn_x, mx_x)) + scale_y_continuous(limits = c(mn_y, mx_y)) +
theme_void()
p3 = qplot() + geom_point(aes(x_wols, y_wols, colour = aes_ewols),
size = 1,
alpha = 0.5) + scale_colour_gradientn(
limits = c(0, mx),
trans = 'pseudo_log',
colours = c("blue", "yellow", "red"),
labels = c(' 0', ' 2', ' 5', '10 ', '20 ', '40 ', '60 '),
breaks = c(0, 2, 5, 10, 20, 40, 60)
) + ggtitle('WOLS') + theme_cowplot() +
ylab(combos[2, i]) + xlab(combos[1, i]) + scale_x_continuous(limits = c(mn_x, mx_x)) +
scale_y_continuous(limits = c(mn_y, mx_y)) + theme_void()
p4 = qplot() + geom_point(aes(x_ols, y_ols, colour = aes_eols),
size = 1,
alpha = 0.5) + scale_colour_gradientn(
limits = c(0, mx),
trans = 'pseudo_log',
colours = c("blue", "yellow", "red"),
labels = c(' 0', ' 2', ' 5', '10 ', '20 ', '40 ', '60 '),
breaks = c(0, 2, 5, 10, 20, 40, 60)
) + ggtitle('OLS') + theme_cowplot() +
ylab('') + xlab('') + scale_x_continuous(limits = c(mn_x, mx_x)) + scale_y_continuous(limits = c(mn_y, mx_y)) +
theme_void()
figure = ggarrange(
p1,
p2,
p3,
p4,
ncol = 2,
nrow = 2,
common.legend = T,
legend = 'right'
) + theme(plot.background = element_rect())
print(annotate_figure(
figure,
left = text_grob(combos[2, i], rot = 90),
bottom = text_grob(paste0(combos[1, i], '\n', '\n'))
))
}
Using lines.
map = SOM(gt_t, xdim = 20, ydim = 20)
e_gt = EmbedSOM(gt_t, map)
e_ngd = EmbedSOM(ngd_t, map)
e_ols = EmbedSOM(ols_t, map)
e_wols = EmbedSOM(wols_t, map)
p1 = ggplot() + geom_segment(
aes(
x = e_gt[, 1],
y = e_gt[, 2],
xend = e_ngd[, 1],
yend = e_ngd[, 2]
),
colour = "red",
alpha = 0.3
) + geom_point(aes(x = e_gt[, 1], y = e_gt[, 2]), size = 0.4, alpha = 0.3) +
ggtitle('NGD') + scale_x_discrete(labels = NULL, breaks = NULL) + labs(x = '', y =
'') + scale_y_discrete(labels = NULL, breaks = NULL) + theme_void()
p2 = ggplot() + geom_segment(
aes(
x = e_gt[, 1],
y = e_gt[, 2],
xend = e_wols[, 1],
yend = e_wols[, 2]
),
colour = "red",
alpha = 0.3
) + geom_point(aes(x = e_gt[, 1], y = e_gt[, 2]), size = 0.5, alpha = 0.4) +
ggtitle('WOLS') + scale_x_discrete(labels = NULL, breaks = NULL) + labs(x = '', y =
'') + scale_y_discrete(labels = NULL, breaks = NULL) + theme_void()
p3 = ggplot() + geom_segment(aes(
e_gt[, 1],
y = e_gt[, 2],
xend = e_ols[, 1],
yend = e_ols[, 2]
),
colour = "red",
alpha = 0.3) + geom_point(aes(x = e_gt[, 1], y = e_gt[, 2]), size = 0.5, alpha =
0.4) + ggtitle('OLS') + scale_x_discrete(labels = NULL, breaks = NULL) + labs(x = '', y =
'') + scale_y_discrete(labels = NULL, breaks = NULL) + theme_void()
p1
p2
p3
ggarrange(p1, p2, p3, ncol = 3)
Color points by per-cell mean squared error.
Error = rep(0, nr)
p1 = qplot() + geom_point(aes(e_gt[, 1], e_gt[, 2], colour = Error),
size = 1,
alpha = 0.5) + scale_colour_gradientn(
limits = c(0, mx),
trans = 'pseudo_log',
colours = c("blue", "yellow", "red"),
labels = c(' 0', ' 2', ' 5', '10 ', '20 ', '40 ', '60 '),
breaks = c(0, 2, 5, 10, 20, 40, 60)
) + ggtitle('Ground Truth') + theme_cowplot() +
ylab('') + xlab('') + theme_void()
p2 = qplot() + geom_point(aes(e_ngd[, 1], e_ngd[, 2], colour = unname(ngd_mcr)),
size = 1,
alpha = 0.5) + scale_colour_gradientn(
limits = c(0, mx),
trans = 'pseudo_log',
colours = c("blue", "yellow", "red"),
labels = c(' 0', ' 2', ' 5', '10 ', '20 ', '40 ', '60 '),
breaks = c(0, 2, 5, 10, 20, 40, 60)
) + ggtitle('NGD') + theme_cowplot() +
ylab('') + xlab('') + theme_void()
p3 = qplot() + geom_point(aes(e_wols[, 1], e_wols[, 2], colour = unname(wols_mcr)),
size = 1,
alpha = 0.5) + scale_colour_gradientn(
limits = c(0, mx),
trans = 'pseudo_log',
colours = c("blue", "yellow", "red"),
labels = c(' 0', ' 2', ' 5', '10 ', '20 ', '40 ', '60 '),
breaks = c(0, 2, 5, 10, 20, 40, 60)
) + ggtitle('WOLS') + theme_cowplot() +
ylab(combos[2, i]) + xlab(combos[1, i]) + theme_void()
p4 = qplot() + geom_point(aes(e_ols[, 1], e_ols[, 2], colour = unname(ols_mcr)),
size = 1,
alpha = 0.5) + scale_colour_gradientn(
limits = c(0, mx),
trans = 'pseudo_log',
colours = c("blue", "yellow", "red"),
labels = c(' 0', ' 2', ' 5', '10 ', '20 ', '40 ', '60 '),
breaks = c(0, 2, 5, 10, 20, 40, 60)
) + ggtitle('OLS') + theme_cowplot() +
ylab('') + xlab('') + theme_void()
ggarrange(
p1,
p2,
p3,
p4,
ncol = 2,
nrow = 2,
common.legend = T,
legend = 'right'
)